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Vertical stacking of two-dimensional (2D) crystals, such as graphene and hexagonal boron nitride, 
has recently lead to a new class of materials known as van der Waals heterostructures (vdWHs) with 
unique and highly tunable electronic properties. Ab-initio calculations should in principle provide 
a powerful tool for modeling and guiding the design of vdWHs, but in their traditional, form such 
calculations are only feasible for commensurable structures with a few layers. Here we show that 
the dielectric properties of realistic, incommensurable vdWHs comprising hundreds of layers can 
be calculated with ab-initio accuracy using a multi-scale approach where the dielectric functions 
of the individual layers (the dielectric building blocks) are coupled simply via their long-range 
Coulomb interaction. We use the method to illustrate the 2D-3D dielectric transition in multi-layer 
M 0 S 2 crystals, the hybridization of quantum plasmons in large graphene/hBN heterostructures, 
and to demonstrate the intricate effect of substrate screening on the non-Rydberg exciton series in 
supported WS 2 . 


The class of 2D materials which started with graphene 
is rapidly expanding and now includes metallic and semi¬ 
conducting transition metal dichalcogenidesQ] in addi¬ 
tion to group III-V semi-metals, semiconductors and 
insulators [2]. These atomically thin materials exhibit 
unique opto-electronic properties with high technologi¬ 
cal potential[3H7]. However, the 2D materials only form 
the basis of a new and much larger class of materials 
consisting of vertically stacked 2D crystals held together 
by weak van der Waals forces. In contrast to conven¬ 
tional heterostructures which require complex and ex¬ 
pensive crystal-growth techniques to epitaxially grow the 
single-crystalline semiconductor layers, vdWHs can be 
stacked in ambient conditions with no requirements of 
lattice matching. The latter implies a weaker constraint, 
if any, on the choice of materials that can be combined 
into vdWHs. 

The weak inter-layer binding suggests that the indi¬ 
vidual layers of a vdWH largely preserve their original 
2D properties modified only by the long range Coulomb 
interaction with the surrounding layers. Turning this 
argument around, it should be possible to predict the 
overall properties of a vdWH from the properties of the 
individual layers. In this Letter we show that this can 
indeed be achieved for the dielectric properties. Concep¬ 
tually, this extends the Lego brick picture used by Geirn 
and Grigorieva [8] for the atomic structure of a vdWH, 
to its dielectric properties. Specifically, we develop a 
semi-classical model which takes as input the dielectric 
functions of the individual isolated layers computed fully 
quantum mechanically and condensed into the simplest 
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possible representation, and couple them together via the 
Coulomb interaction, see Figure[lJ Despite the complete 
neglect of interlayer hybridization, the model provides 
an excellent account of both the spatial and dynamical 
dielectric properties of vdWHs. The condensed represen¬ 
tation of the dielectric function of all isolated 2D crystals 
can thus be regarded as the dielectric genome of vdWHs. 

In addition to its conceptual value, our approach 
overcomes a practical limitation of conventional first- 
principles methods. Such methods are not only compu¬ 
tationally demanding, but also rely on periodic bound¬ 
ary conditions which are incompatible with the incom¬ 
mensurable interfaces found in vdWHs. In fact, for 
many purposes, an in-plane lattice mismatch between 
neighbouring 2D crystals is preferred because it reduces 
the interlayer coupling and thus minimises the risk of 
commensurate-incommensurate transitions [9], and for¬ 
mation of Moire patterns (TO) and associated band struc¬ 
ture reconstructions (TT| which are typical for systems 
with similar lattice constants. This emphasises the need 
for alternative approaches for modelling vdWHs. 

The dielectric function is one of the most important 
material response functions. It determines the effective 
interaction between charged particles in the material, 
contains information about the collective oscillations of 
the electron gas (plasmons) [12], and enters as a funda¬ 
mental ingredient in many-body calculations of e.g. ex- 
citons and quasiparticle band structures 131 ITi]. 

The (inverse) dielectric function is related to the elec¬ 
tron density response function, %, via 

e _1 (r, r ',w) = <5(r-r') + J X {r", r', u)dr”. (1) 

In our quantum-electrostatic heterostructure (QEH) 
model the calculation of the dielectric function is divided 
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FIG. 1: Schematic of the QEH model, (a) The density response function and dielectric function of the 
heterostructure are calculated from the dielectric building blocks of the individual layers assuming a purely 
electrostatic interaction between the layers. The dielectric building blocks are calculated ab-initio for the isolated 
layers. They comprise the monopole and dipole components of the density response function, Xm/Di together with 
the spatial shape of the electron density, Pm/d(z ), induced by a constant and linear applied potential, respectively, 
(b) Monopole and dipole induced densities (blue) together with the associated potentials (red) for monolayer M 0 S 2 . 


into two parts. In the first part the in-plane averaged 
density response function of each of the freestanding lay¬ 
ers, Xii z > z 'i Q||, w ), are obtained from ab-initio calcula¬ 
tions. In practice we treat the in-plane momentum trans¬ 
fer, qy, as a scalar since most 2D materials are isotropic 
within the plane. From \i we calculate the magnitude of 
the monopole/dipole component of the density induced 
by a potential with a constant/linear variation across the 
layer and in-plane variation expiry ■ qy): 

X*a(q||,w) = y z a Xi{z,z'i<i h uj)z ,a dzdz'. (2) 

Here a = 0,1 for the monopole and dipole components, 
respectively. In addition we calculate the spatial form of 
the induced density, Pi a {z, qy). With a proper normal¬ 
ization of pi a we can then write 

J Xi{z,z',<\_\\,u)z lol dz' = Xio(q||,w)p ia ( 2 ,qy) (3) 

We have found that while Xia depends strongly on fre¬ 
quency, pi a does not. The data set ( XionPia .) with 
a = 0,1 or equivalently a = M,D constitutes the dielec¬ 
tric building block of layer i, as illustrated in Figure |Tj 
According to Eq. ([ 3 ]) the dielectric building block allows 
us to obtain the density induced in the (isolated) layer 
s by a constant/lincar potential. It is straightforward 
to extend the dielectric building blocks to account for 
higher-order moments in the induced density described 
by a > 1, but we have found the dipole approximation 
to be sufficient in all cases considered. 


In the second part of the QEH model, the den¬ 
sity response function of the vdWH in the discrete 
monopole/dipole representation is obtained by solving 
a Dyson equation that couples the dielectric building 
blocks together via the Coulomb interaction. The Dyson 
equation for the full density response function giving the 
magnitude of the monopole/dipole density on layer i in¬ 
duced by a constant/linear potential applied to layer j. 
reads (omitting the qy and ui variables for simplicity) 

Xiap/3 = Xia^iaJ0 T Xia E Xkj,jp‘ (4) 

k^i, 7 

The Coulomb matrices are defined as 

^ia,kj (Q|| ) J Pia (z ; Qj|) ^k^ (z > Qjj )dz (5) 

where X.-y is the potential associated with the induced 
density, pfc 7 , which we calculate on a uniform grid by 
solving a ID Poisson equation. Note that we leave out 
the self-interaction terms in Eq. ([ 4 ]) since the intralayer 
Coulomb interaction is already accounted for by the un¬ 
coupled Xia■ The (inverse) dielectric function of Eq. ([I]) 
in the monopole/dipole basis becomes 

W ) = fiia,jf) + ^ (q||) Xfc7,jff(q||, kj). 

k-y 

( 6 ) 

More details on the method and computations are pro¬ 
vided in the supporting information. 
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FIG. 2: The static dielectric function e(qy,w = 0) of the 
51 transition metal dichalcogenides and oxides included 
in the database, that span a large range of magnitudes. 
The relation to the quasiparticle GoWo band gap of the 
materials, calculated in Ref. m is shown in color. 


A database containing the dielectric building blocks of 
a large collection of 2D materials has been constructed, 
and is available from our website m ■ It presently con¬ 
tains more than 50 transition metal dichalcogenides and 
oxides, graphene at different doping levels, and hBN, and 
more materials are being added. From here the data 
files can be downloaded together with a Python mod¬ 
ule for calculating the dielectric function and associated 
properties of any combination of these materials. QEH 
model calculations for vdWHs containing a few hundred 
layers can be performed on a standard PC. In Figure 
[2 we show the qy-dependent static dielectric functions 
of the monolayer transition metal dichalcogenides and 
-oxides presently contained in our database (for a com¬ 
plete overview of the materials see Ref. [16]). All the 
dielectric functions show the same qualitative form, in 
particular they become 1 for qy —>• 0 and qy —> oo, how¬ 
ever there is quite some variation in their magnitude. As 
expected the size of the dielectric function correlates well 
with the size of the band gap of the material indicated 
by the colour. 

First-principles calculations were performed with the 
GPAW code[13, ITS] , Single-particle wave functions and 
energies were calculated within the local density approx¬ 
imation (LDA) using 400 eV plane wave cut-off and at 
least 45 x 45 sampling of the 2D Brillouin zone. Den¬ 
sity response functions and dielectric functions were cal¬ 
culated within the random phase approximation (RPA). 
The RPA does not include electron-hole interactions, but 
generally yields good results for the static dielectric prop¬ 
erties of semi-conductors and dynamical response of met¬ 
als. Except for M 0 S 2 bulk, we included at least 15 A of 
vacuum in the super cells perpendicular to the layers and 
applied a truncated Coulomb kernel to avoid long range 


screening between periodically repeated structures. All 
response functions were calculated in a plane wave basis 
including reciprocal lattice vectors up to at least 50 eV. 
A similar cut off was used for the sum over empty states 
and convergence was carefully checked. The frequency 
dependence of the response functions was represented on 
a non-linear frequency grid ranging from 0 to 35 eV, with 
an initial grid spacing of 0.02 eV. All details of the cal¬ 
culations and atomic structure geometries are provided 
in the supporting information. 

As a first application of the QEH model, we study how 
the (static) dielectric function of a 2D material evolves 
as the layer thickness increases towards the bulk. One 
of the most characteristic differences between 2D and 3D 
materials is the behaviour of the dielectric function in 
the long wave length limit: For a bulk semiconductor, 
the dielectric function e(qy) tends smoothly to a constant 
value larger than unity as qy —► 0. In contrast e(qy) = 
1 + O(qy) for a 2D semiconductor implying a complete 
absence of screening in the long wave length limit mm- 

Ab initio calculations were performed for the dielec¬ 
tric function of M 0 S 2 monolayer, bilayer, and bulk, and 
the QEH model was used for multilayer structures up to 
100 layers. Figure [3] (b) shows the dielectric functions 
averaged over the slabs, i.e. the macroscopic dielectric 
function, as function of the in-plane momentum trans¬ 
fer. For large qy the dielectric functions show similar 
behavior. However, whereas e(0) = 14 for the bulk, the 
dielectric functions of the slabs decrease sharply to 1 for 
small qy. This demonstrates that the dielectric proper¬ 
ties of a vdWH of thickness L are 2D like for qy -C 1 /L 
and 3D like for qy 1 /L. Interestingly, also the result 
for bulk M 0 S 2 shows reminiscence of the 2D nature of the 
constituent layers, where the magnitude of the dielectric 
function has a slight drop when qy —> 0. 

The QEH model describes the change in the dielectric 
function from mono- to bilayer very accurately in spite 
of the well known differences between the mono- and bi¬ 
layer band structures|21|. This shows that hybridisation 
driven band structure effects, i.e. quantum confinement, 
have negligible influence on the dielectric properties of a 
vdWH and is the main reason for the success of the QEH 
model. The model result seems to converge towards the 
ab initio bulk result, however, convergence is not fully 
reached even for N = 100. The slow convergence to¬ 
wards the bulk result is mainly due to the spatial varia¬ 
tion of the induced potential across the slab. In Figure[3] 
(c) we show the z-dependent dielectric function defined 
as e(z) = V e xt/V to t(z) , for a constant (along z) exter¬ 
nal potential with a long wavelength in-plane variation 
for N = 50. Although e(z) is close to the ab-initio bulk 
value (dashed line) in the middle of the slab, screening is 
strongly suppressed in the surface region. Increasing the 
slab thickness beyond 50 layers brings the QEH result 
even closer to the bulk result in the middle of the slab, 
but a small underestimation remains originating from the 
difference in the band structures of the monolayer and 
bulk systems. The suppressed screening in the surface 
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FIG. 3: 2D-3D transition of the dielectric function, (a) Atomic structure of a 20 layer M 0 S 2 slab, (b) The 
macroscopic static dielectric function cm( qy,w = 0) as a function of the in-plane momentum transfer for different 
number of layers, N. The macroscopic dielectric function relates the total potential averaged over the width of the 
slab to an external potential of the form V ext (r\\, z) = exp(«r'|| q||) - The dielectric functions increase monotonically 
with N converging slowly towards the dielectric function of bulk M 0 S 2 obtained from an ab-initio calculation. 
Excellent agreement between the QEH model and the ab initio results are seen for TV = 1, 2. The slow convergence 
towards the bulk result is due to the strong spatial variation of the induced potential in the surface region of the 
slabs. This can be seen in panel (c) which shows V ext /V tot (z), i.e. the local dielectric function, for an external 
potential constant across the slab and with in-plane wave vector qy = 0.036A -1 for N = 50. 


region is a direct consequence of the anisotropic nature 
of the layered M 0 S 2 crystals which limits the screening 
of perpendicular fields relative to in-plane fields, and is 
expected to be a general property of vdWHs. 

The model can also be used to calculate the response to 
fields polarized along the ^-direction, i.e. perpendicular 
to the layers. In this case the perpendicular component, 
e zz (q z = 0), can be calculated by applying an external 
potential with a linear variation along z. In the discrete 
basis of the QEH model, such a field is represented by 
a vector with 0 for all monopole components and 1 for 
all dipole components. Comparing the averaged slope 
of the total potential to the slope of the applied linear 
potential for a slab of N=100 layers of MoS2 yields e 22 = 
7.8. This value is somewhat larger than the bulk value of 
6.03, however, due to long range surface effects it is not 
necessarily to be expected that the two numbers should 
coincide. In fact, we find excellent agreement between 
the QEH model and full ab-initio calculation of e 22 for a 
four layer M 0 S 2 slab (see supporting information). 

Next, we consider the hybridisation of plasmons in 
graphene sheets separated by hBN, see Figure[4]]a). Plas¬ 
mons in graphene on hBN were recently found to prop¬ 
agate with low loss [5], and the close to perfect lattice 
match between the two layers enables full ab initio cal¬ 
culations for the thinnest heterostructures. Here we use 
doped graphene that has a finite density of states at the 
Fermi level, giving rise to two-dimensional plasmons with 
energies in the regime 0-2 eV. The plasmon energies goes 


to zero in the optical limit, qy —> 0 as is characteristic for 
plasmons in 2D metals [22] [53]. We calculate the effect of 
hBN on the plasmons using the QEH model for up to 20 
layers of hBN and compare to full ab-initio calculations 
for 1-3 layers of hBN. 

To identify the plasmons of the heterostructure we fol¬ 
low Ref. [24]. In brief, we compute the eigenvalues, 
e n (w), of the heterostructure dielectric function for each 
frequency point and identify a plasmon energy, hup, from 
the condition Ree n (wp) = 0, see Figure [4jb). The cor¬ 
responding eigenvector, </> n (uip), represents the potential 
associated with the plasmon oscillation, see panel (c). 
This analysis identifies two plasmons corresponding to 
the symmetric (H—|-) and antisymmetric (H—) combina¬ 
tions of the graphene plasmons as previously found for 
two freestanding graphene sheets [25]. For 1-3 hBN lay¬ 
ers, the QEH model perfectly reproduces the ab-initio re¬ 
sults for the dielectric eigenvalues, plasmon energy, and 
weight. The latter was defined as the area under the 
peaks in the loss function — Ime _1 (qy, w), see panel (b). 
The densities and potentials of the plasmon eigenmodes 
shown in panel (c) are also reproduced fairly accurately 
by the model, where the qualitative differences for the 
induced densities, p(z), are due to the use of a limited 
basis of the monopole and dipole response for each layer. 
In panels (e-f) the result of full ab-initio calculations are 
shown by symbols while the QEH results are shown by 
continuous lines. The effect of the hBN buffer (dashed 
lines) is to red shift and damp the plasmons compared to 
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FIG. 4: Plasmons in graphene/hBN heterostructures, (a) Two graphene sheets separated by three layers of hBN. (b) 
Eigenvalues of the heterostructure dielectric function e(u>). Only the two eigenvalue curves that fullfill the plasmon 
condition Ree n (ayp) = 0 are shown, (c) The eigen-potential, i^(wp), and associated density, p(ujp ), of the plasmon 
modes. The plasmons correspond to the antisymmetric (H—) and symmetric (++) combinations of the isolated 
graphene plasmons. (d) Plasmon dispersion for heterostructures containing 1 and 3 layers of hBN. Full lines denote 
the QEH model while ab-initio results are denoted by symbols, (e+f) Energy and weight of the plasmon modes for 
up to 20 layers hBN between the graphene sheets. Results for equivalent structures with vacuum filling the gap are 
also shown. Dashed black lines indicate the plasmon energy and weight in an isolated graphene sheet. Overall, the 
QEH model is in excellent agreement with the full ab initio calculations performed for up to 3 layers hBN. 


the result for two graphene sheets separated by the same 
amount of vacuum (full lines). This is also reflected by 
the relatively large amount of electron density located on 
the hBN during the plasma oscillation, see panel (c). 


Finally, we explore some characteristic features of ex- 
citons in freestanding and supported 2D semiconductors. 
A straight forward generalisation of the well known Mott- 
Wannier model[3?| leads to the following eigenvalue equa¬ 


tion for the excitons of a 2D semiconductor [T51 [28] : 


V 2 

v 2 D 

2/^ex 


W{y) 


F(r)=E b F(r ), 


(7) 


where E b is the exciton binding energy, F( r) is the wave 
function, p ex is the effective mass, and W{ r) is the 
screened electron-hole interaction. Assuming that the 
electron and hole are localised in layer 1, the Fourier 
transformed screened electron-hole interaction is ob¬ 
tained from the static (w = 0) response function Eq. 
Q and Coulomb interaction matrix Eq. ([5| of the QEH 
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FIG. 5: Excitons in supported WS 2 . (a) Monolayer WS 2 adsorbed on a h-BN substrate, (b) The screened 
interaction between an electron and a hole localised within a WS 2 monolayer adsorbed on hBN. For comparison the 
unscreened 1/r potential is shown. The radial probability distribution of the first five excitons, r|F(r)| 2 , are also 
shown (arbitrary normalization), (c) The calculated binding energies of the lowest five excitons in freestanding WS 2 
(green dashed) and WS 2 on hBN (blue) and MoS 2 (cyan). Experimental values from Ref. [215] for WS 2 on Si0 2 are 
shown in red. The 2D hydrogen model with a 1/er potential is shown for e = 1.7. (d) The dielectric function of the 
WS 2 layer defined as e(q) = V(q)/W(q), where V(q) and W{q) are the bare and screened interaction in the WS 2 
layer, respectively, (e) The screened interaction in the WS 2 layer as function of log(r). (f) The relative difference 
between the screened interaction in the supported and freestanding WS 2 . Inset shows the relative difference between 

Et for the supported and freestanding WS 2 . 


model, 

W(q||) = Vim, 1M+ ^lMj/3(q||)Xt/3,ia(q||)^a,lM(q||)- 

( 8 ) 

The first term is the bare, i.e. unscreened, electron- 
hole interaction in layer 1 under the assumption that 
the electron and hole densities can be represented by 
the induced monopole density, Pim{z). The second term 
describes the screening from the surrounding layers and 
layer 1 itself. Note that the above equation can be easily 
generalised to describe the screened interaction between 
charges localised in different layers (relevant for indirect 
excitons). 

In Ref. [2(| Chernikov et al. observed a peculiar non- 


hydrogenic Rydberg series for the excitons in a single 
layer of WS 2 adsorbed on a SiC >2 substrate. Here we use 
the QEH model to calculate the screened electron-hole 
interaction within the WS 2 layer from the dielectric func¬ 
tion of the full heterostructure. Since the QEH is applica¬ 
ble only to layered materials we place WS 2 on a 100 layer 
thick slab of hBN which has dielectric constant very sim¬ 
ilar to that of SiC >2 (both around 4). For comparison we 
performed similar calculations using M 0 S 2 as substrate 
(dielectric constant larger than SiC^)- Figure[5|(c) shows 
the five lowest s-excitons calculated from Eq. (m| for both 
freestanding and supported WS 2 . For freestanding WS 2 , 
we obtain E f, = 0.59 eV for the lowest exciton in good 
agreement with previous ab-initio calculations |29j. The 
enhanced screening from the substrate lowers the exci- 
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ton binding energies bringing the entire series closer to 
the experimental values (red), in particular for the hBN 
substrate. 

The dielectric function of the WS 2 layer defined as 
e(q) = V(q)/W(q), where V{q) and W{q) are the bare 
and screened interaction in the WS 2 layer, respectively. 
Figure[5](d) shows that the dielectric function of the sup¬ 
ported WS 2 layer exceeds unity in the qy —> 0 limit. For 
structures of finite width, L, the dielectric function will 
in practice tend to unity for very small qy <C 1/L. Here 
the result have been extrapolated to infinite substrate 
thickness, where e(qy) tends to a value larger than unity. 
This means that the nature of the screening within the 
layer is not strictly 2D because the bulk substrate is able 
to screen the long wave length fields. In real space, the 
screened potentials diverge as log(r) for small r and de¬ 
cay as 1/r for large r, see panel (e). In panel (f) we show 
how the substrate affects W(r): The relative deviation 
from W(r) of the freestanding layer vanishes for small 
and large r but becomes significant at intermediate dis¬ 
tances. As a consequence, the substrate-induced change 
in the exciton binding energy is relatively larger for inter¬ 
mediate exciton sizes. These results clearly demonstrate 
the profound, nonlocal influence of substrates on the di¬ 
electric screening and excitations in 2D materials. 

In conclusion, we have demonstrated that the spatial 
and dynamical dielectric properties of a vdWH can be 
accurately and efficiently obtained from the dielectric 
properties of its constituent 2D crystals. The presented 


quantum-electrostatic heterostructure model (QEH) ex¬ 
ploits this feature and enables the calculation of the di¬ 
electric properties and collective electronic excitations 
of realistic incommensurable heterostructures with ab- 
initio precision. The dielectric building blocks for more 
than fifty different 2D materials are available in an open 
database allowing 2D materials researchers to efficiently 
predict and design the dielectric properties of realistic 
vdWHs. 
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In this supplementary material we provide a detailed 
description of our quantum-electrostatic heterostructure 
(QEH) model including the precise definition of the di¬ 
electric building blocks. In addition we detail the spectral 
analysis used to identify the plasmon eigen modes for the 
graphene/hBN structures and describe the calculation of 
the screened electron-hole interaction used in the 2D ex- 
citon model. Finally, we provide computational details 
for all the ab-initio calculations presented in the Letter. 


I. FORMAL MATTERS 

Within linear response theory, the induced density due 
to an external field of the form V ext (v, t ) = V ext (r, io)e lut , 
is described by the density response function, x(r,r',w): 



The density response function can be obtained from its 
non-interacting counterpart, x°( r ! r, ,w), that gives the 
response to the total field, by solving the Dyson equation 
in the random phase approximation (RPA): 

X(r,r',w) = x°(r,r',w) + 

dr 1 dr 2 X°(r,r 1 ,u})- --- T x(r 2 ,r , ,w). (2) 

Ti - r 2 

For modelling of vdWHs, this equation is favourably split 
into two parts, namely the intra-layer and inter-layer 
parts, as described below. 

We are assuming a basis set consisting of layer centred 
functions, {4>i a }, where i denotes the layer. Defining the 
Coulomb matrix as 

Vinjfl = J dvdr' 4> ia (y) — 4> jl3 ( r') (3) 

we can divide the Coulomb interaction into its intra- and 
interlayer parts: V = V+V 1 . The Dyson equation 2 can 
then be separated into the following two matrix equations 

X = X° + X°Vx (4) 

X = X + xV J X- (5) 


To see this, simply insert Eq. 4 into Eq. 5 

X = X° + X°Vx + X°V I x + ^VxV 1 * (6) 

= X° + X°V I x + X°V(x + xV I x) (7) 

= X 0 + X°V I x + X°Vx (8) 

= X° + X°(V + V I )X, (9) 

which is the original Dyson equation. 

At this point no approximations, except for the RPA, 
have been introduced. In particular, x° i n Eq. 4 is the 
non-interacting response function of the full vdWH. To 
make progress we make the assumption that the over¬ 
lap/hybridization between wave functions (not to be con¬ 
fused with the basis functions) on neighbouring layers 
can be neglected. This allows us to replace x° of the het¬ 
erostructure by the sum of x° for the individual isolated 
layers. In practice this means that Eq. 4 can be solved 
for each layer separately. 

We calculate x° for the isolated layers within the RPA 
using single-particle wave functions and energies from 
density functional theory (DFT) as described in Ref. 1 . 
The interacting density response function, x, for the 
monolayer is obtained by solving the Dyson equation in 
a plane-wave basis with a 2D truncated Coulomb Kernel, 

Vg'- 

~ 4-7T 

V* d Gz = ^[1-cos(G z L/2)\. (10) 

The use of a truncated Coulomb interaction is essential to 
avoid interaction between periodically repeated layers 2 . 
The truncation length is set to half the unit cell height, 
L. In the plane wave basis, the Dyson equation for the 
density response function, x, is then written: 

XG,G'(qn,w) = XG,G'(qn. w ) + 

X! Xg.Gj (q|| , w) Vgf (q|| )XG 1 ,G' (q|| , w), (11) 

G! 

where qy belongs to the 2D Brillouin zone. 

II. QEH MODEL 

A. The dielectric building blocks 

We start by defining the density response function for 
the individual layers, where the macroscopic average is 




2 



FIG. 1: Basis functions used to represent potentials 
(left) and induced densities (right) in the QEH model. 
The example is for graphene at qj = 0.029A 


taken in the parallel directions. The response function is 
then expressed in terms of the perpendicular coordinates 
^ and z ', and the magnitude of the momentum transfer 
parallel to the layer, gy (we assume isotropic materials, 
where the response does not depend on the direction of 
5 ||, but the method can be straightforwardly generalized 
to non-isotropic 2D materials): 


X{z,z',q\\,uj) 



dr H dr ||X(ry,g||,w) 


1 

L 


G z ,G' Z 


( 12 ) 


and in calculating the matrix elements of the intralayer 
response function we integrate over all space: 



dzdz'<j> i>a (z)x(z,z',qi 


u)<t>i, a {z') 

(16) 


-// dzdz'(z-Zi) a Xi(z,z',q\\,u)(z' -Zi) a , (17) 

where a = {M, D} or equivalently a = {0,1}. 

The basis functions can be interpreted as potentials 
that act on \. In order to represent the induced densities 
produced by these potentials, we introduce two density 
basis functions defined as 


PiA z ’<l\\) 


f dz'x{z, z', g\\,u = 0) 4>i, a ( z ') 

X*,a(?||,w = 0) 


(18) 


As an example, the monopole and dipole density ba¬ 
sis functions for monolayer graphene are shown in 
Fig. 1 (right). We have found that the frequency depen¬ 
dence of the basis functions can in general be omitted, 
while the q\\ -dependence is not always negligible. Divid¬ 
ing by Xi,a{q ||,w = 0) in Eq. 18 ensures that the density 
basis function is normalized such that the overlap with 
the potential basis is unity: (0i, a |p?:,a(i7||)) = 1, where in¬ 
tegration over z is implied. To ease the derivation of the 
Dyson equation in the monopole/dipole basis, we make 
the approximation that the potential and density basis 
functions form a dual basis, i.e. 


where the integration is over the in-plane coordinates, A 
is the in-plane area of the supercell, and L is the height 
of the supercell perpendicular to the layer. Integrating 
over the in-plane coordinates corresponds to taking the 
zero components Gy = Gy = 0 in the plane-wave repre¬ 
sentation of xg,G' (<!:<*>). Working with % instead of \° 
ensures that local held effects within the isolated layer 
are exactly taken into account. 

For an efficient representation of the response functions 
and solution of the Dyson equation we need a small yet 
accurate basis set to represent the induced densities in 
the layers and the potentials created by these induced 
densities. To represent potentials we simply use a con¬ 
stant and linear potential corresponding to a first order 
expansion of the induced potentials, see Fig 1 (left). We 
refer to these as monopole (M) and dipole (D) potentials. 
The potential basis functions of layer i at position Zi are 
thus 


where a,/3 = {M,D}, and i,j are layer indices. This 
implies that, within the subspace spanned by the basis 
functions, we have the completeness relation 

P = El^>«)(0i,a| =1 (20) 

i,OL 

We note that Eq. (19) is not exact because of the small 
but finite overlap between potential and density basis 
functions on neighbouring layers. However, taking this 
into account gives very small modifications to the re¬ 
sulting vdWH dielectric properties. Finally, we note 
that working with a dual basis is natural as, in gen¬ 
eral, the spectral representation of the dielectric func¬ 
tion is written in a dual basis of potential and density 
eigenfunctions 3 . 


II 

TT 

X 

1 [zi — d/2,Zi-\-d/2] 

(13) 

4>i,D ( z ) = 

( z — Zi) 1 [z{ — d/2,Zi~hd/2] 

(14) 

f 1 

if Z 6 C 

(15) 

lc Ho 

if Z $ C 


where d is a localisation parameter that is set equal to the 
interplane distance. Since the density response is already 
confined to the layer, the precise value of d is not essential 


B. Electrostatic Dyson equation 

The Dyson equation (5) for the heterostructure den¬ 
sity response function x( z i z ’ > Q \\, w) is now written in the 
potential basis of dimension 2 N x 2 N, where N is the 
number of layers. In the following the (gy, w) variables 
are omitted from the expressions for simplicity. Response 
functions x, \ and Coulomb kernel V are regarded as op¬ 
erators and integration over r, r' is implied in the inner 
























3 


products. The matrix elements of x are written in the 
potential basis: 

(<Pi,oc]x\<Pj,d) = ( 0 i,a|xl 0 j„ 3 > + ( 0 i.a|x^ 7 xl^,/ 3 >- ( 21 ) 

The first term on the right hand side is simply the re¬ 
sponse function of the isolated layers for which we have 
(‘kalxl^/s) = Xi, a 5 ia ,j0- In the second term, applying 
(0i,a| to x returns Xi,a(Pi,a| (this follows from Eq. 18 
and the symmetry of x(z, z')). Now the completeness 
relation (20) is inserted between V 1 and x, leading to 

(^i.alxl^j,/?) = Xi,a^iaJ0 T 

Xi.c'EiPijV 1 |Pfc,Q! / )( 0 fc,Q: / Ix|0j,j9) 

k,Oi' 

This leads to the final Dyson equation for the het¬ 
erostructure: 

i.Q\\ 5 Xi,ct {,Q\\ ? (3 T 

Xi,a.(Q || j^) ^ ] Via,ky{Q\\) Xk-y,j/3(Q\\ ? ^)- (22) 

k^i ,'y 

The Coulomb kernel is here defined in the density basis 
as: V)a,fca' = {pi,a\V\pk,a')- The term V\Pk, a ') is the 
potential at z from the density basis function in layer k, 
which is found by solving Poisson’s equation for | pk, a ') 
on a real space grid. Since the density parallel to the 
layer just shows periodic oscillations with wave vector gy, 
Poisson’s equation reduces to a ID differential equation: 

d 2 

(-^0 tf|| *^kot' (z) — 47 Tp} ca '(z'). (23) 

The elements of the V matrix are then: Vi a ^ a ' — 

{Pi,a | *&k,a ') ■ 

C. The dielectric matrix 

The inverse dielectric function is related to x through: 
c _1 = 1 — Vx- Due to the non-symmetric nature (in r 
and r') of the dielectric function, the elements of e -1 are 
naturally written using a mixed density/potential basis: 

(p?:,a| e Ifij, a) = fiictJ/3 T (pj,a|R Xl^j,^)- (24) 

Upon insertion of the completeness relation (20) this 
gives 

e ia,jt 3(911 1 W ) = ^ ' Via,k~/{< 1 \\ )Xk'y,jfi{q\\ ; &)■ (25) 

fc,-y 


III. PLASMONS EIGENMODES 

By following a previously developed method for iden¬ 
tifying plasmon eigenmocles in nanostructures from ab 


initio 3 , the dielectric matrix for the heterostructure, 
Eq. 25, is diagonalized to solve the eigenvalue equation: 

^ 2 £ioc,j 0 (q\\,v)fn,j 0 (q\\,u) = e n (q\\,w)f ntia (q\\,u)), 
id 

(26) 

which returns the eigenvalues, e„(gy,w), and eigen¬ 
vectors, fn,ia(q\\,w) °f the dielectric matrix in the 
monopole/dipole basis. A plasmon eigenmode fullfills 
that: 

Re e ia,j 0 (q \\, uJ )f n , j 0 (q \\, w ) = o, (27) 
id 

corresponding to Ree„(gy,a;) = 0. In practice, the plas¬ 
mon energies are identified from the peaks in the eigen¬ 
value loss-spectrum —Ime n (gy,a;) since this includes the 
finite imaginary part which can shift the plasmon energy. 
The right eigenfunctions f n ,ia give the induced potential 
of the plasmon in the basis of (j>i t M/D■ The left eigen¬ 
functions, f '" a , correspond to the induced density of the 
plasmon in the basis of P^m/d 3 - The induced density is 
thus given by 

Pn{z,Q\\) -’52,f? a Pi<x{z,q\\) (28) 

io 

with the corresponding induced potential 

<t>n(z,q\\) ~ yi fr a ®ia(z,q\\) (29) 


IV. EXCITONS 

The Mott-Wannier model, widely used to model ex- 
citons in bulk semiconductors, can be straightforwardly 
generalised to 2D semiconductors. This leads to a 2D 
hydrogenic Hamiltonian of the form 

-p^+W(r)]F(r) = E b F(r), (30) 

Z/lex 

where F{ r) is the exciton wave-function, 0 ex the exciton 
effective mass and W(ry) is the screened Coulomb po¬ 
tential which includes the screening coming from the 2D 
material itself and the environment, e.g. a substrate. 

Now, consider electron and hole charge distributions 
given by (the in-plane variation is a plane wave of wave 
vector gy) 

p e/ T>9u) = YjP^ h ( q \0 pia ( z ' q \0 ( 31 ) 

iot 

We can then calculate the screened interaction between 
the electron and hole charge distributions according to 

W(q\\)= E pka(q\\) e ka,id( , i\\)Vid,ji(.q\\)Pj 1 (q\\)- 

(32) 
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In the case of excitons located in the layer 1 we can 
approximate p e ' h (z,q\\) = piM(z,q\\)(z) and we recover 
the expression in the Methods section. We can de¬ 
scribe a general charge distribution, e.g. using con¬ 
duction/valence band charge distributions p e ' h (z,q\\) = 
W’c/viz, q \\)\ 2 , by a simple redefinition of the Coulomb 
matrix elements in Eq. 32. 

Performing a 2D Fourier transform of W (g||) yields the 
screened potential in real space: 

w ( r \\) = ~^ K J Q d 1 \m J o(Q\\ r \\)W{q\\), (33) 

where J„(x) is the Bessel function of the first kind. The 
exciton mass can be obtained e.g. from an ab-initio band 
structure calculation. We solve Eq. 30 using polar coor¬ 
dinates and a logarithmic radial grid. 


V. SCREENING OF PERPENDICULAR FIELDS 

In Fig. 3c in the manuscript, the spatial form of the 
response due to a constant perturbation across a N=50 
layer slab of M 0 S 2 is shown. This gives the dielectric 
function due to a finite wavevector in the plane of the ma¬ 
terial. However, the model can also be used to calculate 
the response to fields with a variation in the z—direction, 
perpendicular to the layers, and can thus be used to cal¬ 
culate the z—component of the dielectric function, e zz . 
This can be calculated in the optical limit, q z —> 0 with 
the expression: 

L/2 

—L/2 

where L is the width of the structure. In the heterostruc¬ 
ture model, this corresponds to taking the matrix prod¬ 
uct of e~ajp{q\\, w) with a vector, v, with the elements: 
Vjp = where only the dipole elements are non-zero: 
v = {0,1, 0,1,...}. The expression becomes: 

fQ 1 = jj Vi ’ D e iD,jD(<l\\ = °, u = °) v i,D, (35) 
ij 


//: 


ze : (z, z'oj = 0) z' dz'dz, (34) 


< 


o 

a 


H 


Perpendicular response, M 0 S 2 



FIG. 2: Induced potential of a N=4 layer M 0 S 2 slab, 
due to an external perturbation with a constant slope 
across the structure. The potentials are normalized 
with respect to the potential drop across the structure, 
AI4 xt , the width of the structure here defined as 
L = 4dMoS 2 = 24.6A. The dashed lines indicate the 
center of the outermost layers. 


that determines the total potential. The induced poten¬ 
tial is then obtained as Vj n d = Vtot — I4 xt . 

As seen in Fig. 2, the QEH model captures the re¬ 
sponse to perpendicular fields quite well, with a ten¬ 
dency to overestimate the drop in induced potential 
across the structure and therefore overestimate the di¬ 
electric function. This leads to a value of e zz (QEH) = 
7.71 compared to an ab initio value of e zz (ab initio) = 
6.81 for the N = 4 M 0 S 2 slab. In case of bulk 
M 0 S 2 we obtain e 22 (bulk, ab initio) = 6.03 compared to 
e zz (N = 100, QEH) = 7.83, which means that the bulk 
limit is less well-described. However, this is to be ex¬ 
pected since the model cannot account for the bulk limit 
as g|j —> 0, since the dielectric function e(<jj| —> 0) = 1 for 
finite slab widths in the model, while for a 3D system the 
dielectric function tends to a finite value. 


VI. COMPUTATIONAL DETAILS 


where N is the number of layers. 

I 11 Fig. 2 the induced potential of a N=4 layer M 0 S 2 
slab due to an external potential with a linear form along 
z, I4 x t(z) oc z is shown together with the ab initio result. 
The potential is clearly screened by the material, where 
the induced potential has opposite sign that the external 
potential. The ab initio result is in this case obtained by 
applying a weak electric held (within the linear response 
regime) in the z-direction on the ground-state DFT level. 
This calculation was performed on a real-space grid rep¬ 
resentation of the electronic wavefunctions, with a grid¬ 
spacing of h = 0.18 A, and (12,12) fc-points, which were 
sufficient to converge the ground-state electronic density 


A. Multilayer M0S2 

Ab initio calculations were performed for monolayer 
M 0 S 2 to obtain the monolayer density response functions 
and induced densities used as input for the heterostruc¬ 
ture model. The single-particle energies and wave func¬ 
tions were calculated with the PBE exchange correlation 
functional, with a plane-wave basis set with an energy 
cutoff of 400 eV. A dense fc-point sampling of (128,128) 
in the 2D Brillouin zone was used in order to calculate 
the response at low momentum transfers. In the linear 
response RPA calculation we used an energy cutoff of 50 
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eV for the reciprocal lattice vectors. We used a nonlinear 
frequency grid from 0 to 35 eV, with an initial grid spac¬ 
ing of 0.02 eV and a broadening of 0.04 eV. Correspond¬ 
ing ab initio calculations were performed for bulk and 
bilayer M 0 S 2 , but with a fc-point sampling of (64,64,1) 
for the bilayer and (64,64, 8) for bulk. For the monolayer 
and bilayer calculations the truncated Coulomb kernel, 
see Eq. 10, was used while the full, i.e. non-tmncated 
kernel, was used for the bulk calculation. We used an 
in-plane lattice constants of 3.18 A, and A-B stacking 
with 6.15 A separation between layers. For the mono- 
layer and bilayer calculation the unit cells contained 20 
A of vacuum to separate the periodic images in the z- 
direction. For the heterostructure calculation we used 
the same separation between the layers as for the ab ini¬ 
tio calculations (d = 6.15 A). We note that the effect of 
stacking arrangement (A-A or A-B) cannot be accounted 
for within the model. 


B. Graphene/hBN heterostructures 

Ab initio calculations were performed to obtain the 
dielectric building blocks of monolayer doped graphene 
and hBN. Also, full ab initio calculations were done for 
entire heterostructures, including up to three layers of 
hBN, or the equivalent amount of vacuum, separating 
the doped graphene layers. An in plane lattice-constant 
of 2.5 A was used for both graphene and hBN, so that 
the heterostructure could be represented a 1 x 1 unit cell. 
The layers were stacked in A-B configuration, with 3.326 
A separation (c-lattice constant of 6.653). We used PBE 
exchange-correlation, a 340 eV energy cutoff for the plane 
waves in the ground state calculations, and (100,100) fc- 
point sampling in the 2D Brillouin zone. In the response 
calculation doped structures were obtained by shifting 
the Fermi-level 1 eV upwards. An energy cutoff of 70 
eV was used for the reciprocal lattice vectors, and un¬ 
occupied bands were included up to 35 eV above the 
Fermi level. All the calculations employed the truncated 
Coulomb interaction and 20 A vacuum to separate the 
repeated structures. A non-linear frequency-grid with an 
initial grid spacing of 0.02 eV and a broadening of 0.05 
eV was used to represent the dynamic response function. 
Plasmon eigenmodes were obtained by diagonalizing the 
dielectric matrix in Bloch representation as described in 
ref. 3 . 


C. Excitons in supported WS 2 

The dielectric building blocks of the WS 2 , hBN, and 
M 0 S 2 monolayers were calculated as follows. Single¬ 


particle energies and wave functions were calculated us¬ 
ing LDA, a plane wave cut-off of 500 eV, and (45,45) 
fc-points. The density response function was calculated 
within RPA using an energy cut-off of 300 eV and includ¬ 
ing empty states up to 50 eV above the Fermi level. The 
truncated Coulomb kernel was employed and 20 A vac¬ 
uum was included in the supercell to separate repeated 
layers. In setting up the heterostructure we used a sep¬ 
aration of 3.22 A between the 100 layers of h-BN and 
5.08 A between WS 2 and h-BN. For WS 2 on 50 layers 
of M 0 S 2 we used a uniform separation of 6.3 A between 
all layers. We then calculated the screened interaction 
from Eq. 32 for gy up to (and including) the second Bril¬ 
louin zone. For calculating the exciton Rydberg series we 
solved Eq. 30 for spherical states on a radial logarithmic 
grid and verified that the exciton energies were converged 
to within 0.01 eV. 


D. 2D Database 

The dielectric building blocks were calculated for 51 
transition metal dichalcogenides and oxides, hBN, and 
graphene at 10 different doping levels from 0.1 to 1 eV. 
For the single particle wave functions and energies ob¬ 
tained from DFT, we used PBE exchange-correlation and 
a plane-wave basis with a energy cutoff equal to 500 eV. 
The 2D Brillouin zone was sampled by (200,200) fc-points 
for graphene, and for the remaining materials we used a 
fc-point density corresponding to (100,100) fc-points. 

For the density response functions we used a cutoff of 
100 eV for the transition metal dichalcogenides and ox¬ 
ides and 150 eV for graphene and hBN. The truncated 
Coulomb kernel was employed and 20 A vacuum was in¬ 
cluded in the supercell to separate the repeated layers. 
All materials were represented on the same frequency grid 
from 0 to 35 eV, with an initial spacing of 0.01 eV and 
a broadening of 0.05 eV. The response functions were 
calculated for a range of in-plane momentum transfers, 
gy, within the first Brillouin zone of graphene up to a 

maximum value of gy = 2.89 A 1 . At small gy below 
0.3 A 1 we use a denser sampling with a grid spacing of 
0.015 A“ in order to capture the strong gy-dependence 
of the plasmon energies and the dielectric function in this 
region. After this limit the grid spacing is increased to 
0.029 A . In order to obtain all response functions on 
the same gy-grid, the data for the remaining materials 
was interpolated to the grid for graphene using conven¬ 
tional 2D spline interpolation. 


* Electronic address: thygesen@fysik.dtu.dk 
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